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Recently, it was argued that a three-point susceptibility equal to the density derivative of the 
intermediate scattering function, Xn{k;t) — dF(k;t)/dn, enters into an expression for the diver- 
gent part of an integrated four-point dynamic density correlation function of a colloidal suspension 
[Berthier et al, J. Chem. Phys. 126, 184503 (2007)]. We show that, within the mode-coupling 
theory, the equation of motion for x n (k;t) is essentially identical as the equation of motion for the 
q — » limit of the three-point susceptibility Xq(k;i) introduced by Biroli et al. [Phys. Rev. Lett. 
97, 195701 (2006)]. We present a numerical solution of the equation of motion for Xn(k\t). We also 
derive and numerically solve an equation of motion for the density derivative of the self-intermediate 
scattering function, Xn(k;t) = dF s (k\t) / dn. We contrast the wave vector dependence of Xn{k;t) 
and x s „ (&;*)■ 

PACS numbers: 64.70.pv, 05.70.Ln 



I. INTRODUCTION 

While it has been generally accepted for quite some 
time that upon approaching the glass transition the liq- 
uid's dynamics become increasingly heterogeneous [l|, 
theoretical descriptions of various four-point correlation 
functions which quantify dynamic heterogeneity have 
been scarce @, Q . It has been only in the last few years 
that a number of different theoretical approaches to this 
problem appeared [jjaja LZI, |8|, |9|, • According to one 
of these approaches [1, [§] , the divergent part of a four- 
point dynamic density correlation function, integrated 
over the whole space, can be expressed in terms of so- 
called three-point susceptibilities. These susceptibilities 
arc defined as the derivatives of the usual two-point corre- 
lation function, i.e. the intermediate scattering function, 
with respect to the densities of the conserved quantities. 
In the case of a colloidal system there is only one con- 
served quantity, the number of particles and thus there 
is only one susceptibility, the density derivative of the 
intermediate scattering function. 

The first goal of this note is to derive the equation of 
motion for this susceptibility within the mode-coupling 
approximation [111 . Il2t Il3j ] and to present numerical solu- 
tions of this equation. Not surprisingly, the equation of 
motion for the the density derivative of the intermediate 
scattering function is essentially identical to that for the 
long-wave vector limit of the three-point susceptibility 
introduced by Biroli et al. [6(. The latter susceptibility 
was defined as the derivative of the intermediate scatter- 
ing function with respect to an external, inhomogeneous 
potential. 

It should be noted that while most theoretical ap- 
proaches consider the full (i.e. collective) four-point dy- 
namic density correlation function, in computer simula- 
tions usually the self part of this function is monitored 
[Tij . For this reason, while discussing the wave vector de- 
pendence of dynamic heterogeneities, Berthier et al. used 
a three-point susceptibility equal to the density deriva- 
tive of the self-intermediate scattering function (see Ref. 



@, Sec. IIC). Parenthetically, in this way they intro- 
duced a somewhat ad-hoc extension of their theoretical 
approach. The approach of Berthier et at, at least as 
presented in Refs. [g, |9j, is only applicable to the full 
four-point function. 

The second goal of this note is to derive and numeri- 
cally solve the equation of motion for the density deriva- 
tive of the self-intermediate scattering function within 
the mode-coupling approximation, and to contrast the 
time and wave vector dependence of the density deriva- 
tives of the collective and self-intermediate scattering 
functions. In Ref. [9( the density derivative of the self- 
intermediate scattering function was calculated by solv- 
ing the mode-coupling equation of motion for this func- 
tion at different densities and then performing a numeri- 
cal differentiation. Numerical solution of the equation of 
motion for the density derivative of the self-intermediate 
scattering function allows us to investigate the wave vec- 
tor at which, according to Berthier et al. [1, the 
dynamics of a glassy colloidal suspension is maximally 
heterogeneous. 



II. EQUATIONS OF MOTION FOR Xn{k;i) AND 

Xn(M) 

We consider the simplest possible model of a colloidal 
suspension: a system of interacting Brownian particles. 
Our first object of interest is a three-point susceptibility 
Xn(k;t) which is defined as the density derivative of the 
intermediate scattering function, 



Xn(k;t) 



dF(k;t) 
dn 



(1) 



Here F(k\ t) denotes the intermediate scattering function, 



= - (n(M)n(-k)) 



(2) 



where N is the number of particles, ri(k; t) is the Fourier 
transform of the microscopic density at time t and n(k) = 



2 



n(k;t = 0). Finally, in Eq. JT]) n is the number den- 
sity, n — N/V, with V being the volume of the system 
(the particle diameter a is used as the unit of length and 
thermodynamic limit is implied throughout). The initial 
value of the intermediate scattering function is equal to 
the static structure factor S(k), F(k; t = 0) = S(k). 

Using standard, exact but formal, methods (e.g. a pro- 
jection operator approach) one can derive the following 
equation of motion for the intermediate scattering func- 
tion, 



dt'(5(t - t') + M irr (fc; t - f))d t ,F(k; t') = 



Dpk 2 
S(k) 



F(k;t). 



(3) 



where Do is the diffusion coefficient of an isolated Brown- 
ian particle and M lrr (fc; t) is the so-called irreducible [l5[ 



memory function. Within the mode-coupling approxi- 
mation [ll|, [l2j applied to the the system of interacting 
Brownian particles, the irreducible memory function is 
given by the following expression [l3| 



M hl (k;t) = 

^/ (§H c(fci)]< - ki 

xF(fci;t)F(|k-ki|;i). 



c(|k-ki|)k- (k-ki; 



(4) 



Here k = k/fc and c(k) is the so-called direct correlation 
function, c(k) = (1 - l/S(k))/n. 

Substituting Eq. J3J into Eq. ((3|) and then differenti- 
ating the resulting equation with respect to the density 
we obtain the following equation of motion for the three- 
point susceptibility x„(fc;i), 



dt'(6(t - + M"(fc;i - t'))d t , X n(k;f) + -^rXn(k;t) 
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S(k) 



-nD 
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D k 2 dS(k) 
S 2 (k) dn 
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F(k;t) 
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\ dn 



dn 



(5) 



We notice that the left-hand-side of Eq. ((5]) is exactly 
the same as the left-hand-side of the overdamped limit of 
the equation of motion derived by Biroli et al. [fj] for their 
three-point susceptibility x q (k;i) in the q — » limit. 
The three-point susceptibility Xq(k;i) is defined as the 
derivative of the intermediate scattering function with 
respect to a static, inhomogeneous external potential. 

Furthermore, we recall that the derivation of the mode- 
coupling expression for the irreducible memory function 



invokes the convolution approximation for the three- 
point equilibrium density correlation function 13]. The 
convolution approximation is equivalent to assuming that 
the direct correlation function is density-independent 
p73 |. To be consistent, we have to use the same ap- 
proximation in Eq. ((5J. Using dc(k)/dn = and 
dS{k)/dn = S 2 {k){c{k)+ndc(k)/dn) = S 2 {k)c{k) in Eq. 
(|5j) we obtain the following equation of motion for the 
three-point susceptibility Xn(k;t) 



" dt'(S(t - t') + M iri (k; t - t'))d t/Xn {k; t') + ^rXn{k; t) 



+nD Q 



S(k) 
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D k 2 c(k)F(k;t) 
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At this point we notice that the difference between the 
right-hand-side of Eq. ^ and the right-hand-side of the 
equation of motion derived by Biroli et al. 6] for the 
three-point susceptibility Xq(k t) in the q —> limit is a 
constant factor equal to —nS(0). This factor is equal to 
the thermodynamic derivative, —nS(0) = — (dn/df3fj,) T 
where (3 — l/(fcsT) and \x is the chemical potential. To 
rationalize this fact we notice that Biroli et al. used ap- 
proximations which in the limit of q — > amount to 
dc(k)/dn = 0. Furthermore, their three-point suscep- 
tibility is defined as the derivative of the intermediate 
scattering function with respect to the external poten- 
tial. In the long wavelength limit, q — > 0, this derivative 
differs from the derivative with respect to the density by 
a thermodynamic factor proportional to [dn/ d^ji) T . 

Our second object of interest is a three-point suscepti- 
bility Xn(k',t) which is defined as the density derivative 
of the self-intermediate scattering function, 



x s n (k;t) = 



dF s {k-t) 
dn 



(7) 



Here F s (k;t) denotes the self-intermediate scattering 
function, 



F s (k;t) = (niCMJmt-k)), 



(8) 



where ni(k; t) is the Fourier transform of the microscopic 
density of one selected (labeled) particle at time t and 
ni(k) = ni(k;t = 0). The initial value of the self- 
intermediate scattering function is equal to 1, F s (k;t — 

0) = 1- 

The derivation of the equation of motion for Xn(k] t) is 
analogous to that for Xn(k;t). Here we only present the 
starting point, i.e. the equation of motion for F s (k;t) 
and the final result, i.e. the equation of motion for 
Xn (k; t). The equation of motion for the self-intermediate 
scattering function F s (k;t) reads 



dt'(6(t - t') + M sirr (fc; t - t'))d t iF s {k; t') 



-D k 2 F s (k;t) 



(9) 



where the self irreducible memory function is given by 
the following expression, 



M silI (k;t) 



(10) 
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The equation of motion for three-point susceptibility 
Xn(k;t) has the following form, 



dt'(S(t - t') + M si "(k; t - t'))d v x s n {k; t') + D k 2 X s n (k; t) 



+nD 
= -D I dt' 
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^3 (c(fci)k • ki) F(h;t - t')F s (\k - ki|; t - t')d t ,F s {k; t') 



(27T) 



(c(fci)k • kij xn(h;t - t')F s (\k - ki|; t - t')d t ,F s {k- 1') 
I 



(11) 



Note that to derive Eq. (fTTj) we again used dc(k)/dn = 0. 



III. NUMERICAL SOLUTIONS OF EQUATIONS 
OF MOTION FOR Xn(k;t) AND x°(fc;t) 

Equations of motion |6]) and (TTTj) for the three-point 
susceptibilities Xn(k;t) and xt (k\ t) can be solved using 
the algorithm used previously jl7l fl8l . [T^ to solve mode- 
coupling equations ((3j4j) for the intermediate scattering 
function. The only input required is the structure factor 
S(k). Here, in order to be consistent with earlier re- 
lated work Q, we use the structure factor calculated for 
the hard sphere interaction potential using the Percus- 
Yevick closure approximation. As is customary for a 
mode-coupling calculation for the hard sphere system, 



we report the results using volume fraction </> = nira 3 /6 
where a is the hard sphere diameter or relative distance 
from the mode-coupling transition e = {<f> c — <j))/(j) c . 

To solve equations of motion ([6|) and (fTT|) we used 300 
equally spaced wave vectors from k — to k — 60 with 
the first wave vector fco = 0.1. For this discretization of 
the mode-coupling equations (0H]) the ergodicity break- 
ing transition (i.e. the mode-coupling transiton) is lo- 
cated at volume fraction <fi c — 0.515866763. We also 
performed a few calculations with larger cutoffs for the 
integral and/or a finer grid of wave vectors and obtained 
qualitatively the same results. 

In Figs. Q]and[2]we show Xn{k;t) and Xn(k',t), respec- 
tively, as functions of both time t and wave vector k for 
two different values of e = (<p c — (f>)/(j> c , e = 0.05 and 
e = 10~ 4 . The former value of e is comparable to the 
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FIG. 1: Three-point susceptibility Xn(k;t) for e = 0.05 (upper 
panel) and e = 10 -4 (lower panel). Contours correspond to 
Xn(k;t) — 4 m where m is an integer, starting from m — — 1. 
The arrow marks the position of the first peak of the static 
structure factor. 



lowest reduced temperature (T — T c )/T c at which mode- 
coupling theory agrees with computer simulations for 
the well-known glass former, the Kob- Andersen Lennard- 
Jones binary mixture [2(J. Using contour plots to show 
these three-point susceptibilities is inspired by a simi- 
lar presentation used by Lechenault et al. [21| in their 
investigation of dynamic heterogeneity in dense, two- 
dimensional granular systems. One should note, how- 
ever, that Lechenault et al. monitored and analyzed a 
correlation function of the fluctuations of the self-overlap 
integrated over the whole space. Here, we present results 
for three-point susceptibilities and we need to invoke the 
theory of Berthier et al. to connect these susceptibilities 
to dynamic heterogeneity. 

We note that, as expected, the susceptibility Xn{k\t) 
has a well-defined maximum at a characteristic time and 
wave vector. The characteristic time strongly increases 
upon approaching the mode-coupling transition whereas 
the characteristic wave vector is approximately indepen- 
dent of the distance from the transition and equal to 
7.1 close to the transition (note that according to the 
Percus-Yevick approximation the first peak of the static 
structure factor at 4> c is located at k — 7.1). 

Examining Fig. [5] we note that, somewhat surpris- 
ingly, Xn(k',t) does not have a well-defined maximum at 
a characteristic time and wave vector. In contrast with 
Xn{k',t), x*(A;;i) as a function of time and wave vector 
forms a ridge which, with decreasing wave vector, gently 
increases and moves towards longer times. For a fixed 
wave vector, e.g. for k = fc max = 7.1 (i.e. the position 
of the first peak of the structure factor at </> c ) , the peak 
position of Xn(k]t) strongly increases upon approaching 
the mode-coupling transition. 



FIG. 2: Three-point susceptibility Xn(k',t) for e = 0.05 (upper 
panel) and e = 10 -4 (lower panel). Contours correspond to 
Xn(k; t) — 4" 1 where m is an integer, starting from m — — 1. 
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FIG. 3: Right scale: dependence of the peak time of 
Xn (fcmax; t) and Xn (fcmax; t) , Tmax and r max , respectively, and 
of the a relaxation time [23], r a , on the reduced distance from 
the mode-coupling transition, e = (4> c — 4>)/4> c . Squares: r ma x; 
filled diamonds: r max ; circles: r a . Left scale: dependence of 
the peak value of three-point susceptibilities, x« (fcmax; T max ) 

X n (fcmax , Tmax 

) on e. Open triangles: x« (fcmax ; T ma x); 
filled diamonds: Xn (fcmax; r max ). 



In Fig. [3] we analyze in some detail the dependence 
of both three-point susceptibilities on the reduced dis- 
tance from the mode-coupling transition for a fixed wave 
vector equal to the position of the first peak of the 
static structure at the transition, k = /j max = 7.1. We 
show the dependence on e = (<f> c — <j))/4> c of the time 
at which Xn(k- msx ;t) and Xn(k max ;t) has the maximum 
value, r max and T^ ax , respectively, and we compare these 
times to the a relaxation time [22J, r Q . As expected, 
have the same e dependence as r a , 



both r max and r^ ax 
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FIG. 4: Wave vector dependence of the Xn(^i r max) (solid 
curves) and of the maximum value of Xn(k; t) (dashed curves) 
for e = 0.05 (lower curves) and e — 10 -4 (upper curves). 




FIG. 5: Wave vector dependence of the Xn(fc; r max) (solid 
curves) and of the maximum value of Xn{k; t) (dashed curves) 
for e = 0.05 (lower curves) and e — 10~ 4 (upper curves). 

Tmax ~ T max ~ T a ~ e and ratios of these times 
and T a are approximately constant, T max /r Q w 1.4 and 
T max/ r Q ~ 0.96. In addition, we show in Fig. [3] the de- 
pendence on e of the peak value of the three-point suscep- 
tibilities for a fixed wave vector k = fc max = 7.1. As ex- 
pected, we find x„(A;max;T- max ) ~ X^max; T^ax) ~ e ~ l - 

The original motivation for investigating two- 
dimensional (contour) plots of Xn{k;t) and \n{k',t) was 
to find the characteristic wave vector at which these sus- 
ceptibilities have maximum values. Once these wave vec- 
tors have been found, one could use the theory of Berthier 
et al. to claim that at these wave vectors dynamic het- 
erogeneities are the strongest. From this perspective, 
Figs. [T][2]lead to a rather surprising conclusion: at a fixed 
time comparable to the a relaxation time both Xn(k;t) 
and Xn(k> have a maximum at a nonzero wave vector. 
However, Fig. [2] suggests that the absolute maximum 
of Xn(k',t) is located at or close to k = 0. The differ- 
ence between the wave vector dependence of Xn(k] t) and 
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FIG. 6: Scaling plot of Xn{kmax;t). Solid lines correspond 
to e equal to 10"\ 10" 2 , 10" 3 , 10~ 4 , 10" 5 , and 10" 6 (top to 
bottom). Inset: time dependence of Xn(&max;i) at a very 
small value of e = 10~ 6 . Dashed lines in the main figure and 
in the inset indicate transient power law dependence in the 
early and late (3 regimes. 

Xn(k;t) is investigated in some detail in Figs. [Hand [5] 
In the former figure we compare the wave vector depen- 
dence of Xn(k;T max ) (recall that r max is the peak posi- 
tion of Xn(k max ;t)) with the wave vector dependence of 
the maximum value of Xn{k] t) for two different values of 
e — (cj) c ~ <j>)/(j)c, e = 0.05 and e = 10~ 4 . We notice that 
the wave vector dependence is quite strong but qualita- 
tively similar, except for higher wave vectors. In Fig. [5] 
we show the wave vector dependence of Xn(k', T max ) (r^ ax 
is the peak position of XnC^maxi t)) with the wave vec- 
tor dependence of the maximum value of x«(^; t) for the 
same two different values of e. We see that in this case 
the wave vector dependence is quite a bit weaker than 
that shown in Fig. |U More importantly, we see that the 
wave vector dependence of Xn(^i r max) and of the the 
maximum value of Xn(^i^) are qualitatively different. 

Finally, in Fig. [5] we show that upon approaching 
the mode-coupling transition three-point susceptibility 
XnC^max! t) approaches the scaling limit. Three point sus- 
ceptibility Xn{kma.x]t) also approaches the scaling limit, 
but in a somewhat more complicated way (not shown). 
The inset in Fig. [6] shows the time dependence of 
Xn(^max;£) for a small value of e equal to 10~ 6 . Com- 
parison of the main figure and the inset shows that upon 
approaching the mode-coupling transition the time de- 
pendence of both susceptibilities at k — fc max is quite 
similar. 



IV. DISCUSSION 

The main reason for our interest in three-point suscep- 
tibilities Xn(k] t) and Xn(^i was the connection between 
these susceptibilities and four-point correlation functions 
which quantify dynamic heterogeneities. Roughly speak- 
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ing, according to Berthier et al. [1, Q , Brownian system's 
four-point correlation function Xi(k;t), which is defined 
as the self part of the four-point dynamic density correla- 
tion function integrated over the whole space, is propor- 
tional to t) 2 for a system in which the total number 
of particles can fluctuate and is proportional to Xn(k',t) 
for a system with fixed total number of particles (note 
that the latter one is commonly used in numerical sim- 
ulations). Following Berthier et al. and assuming that 
mode-coupling theory gives at least qualitatively correct 
results for three-point susceptibility Xn{k', t) we are faced 
with a striking conclusion that the there is no finite in- 
trinsic |23j wave vector at which dynamic heterogeneities, 
as quantified by x^(k;t), are the largest. More precisely, 
at a fixed time we can determine a finite characteristic 
wave vector, but with increasing time this wave vector 
is decreasing towards 0. In other words, the maximum 
value of Xi(k',t) is predicted to monotonically decrease 
with increasing k. 

The wave vector dependence of four-point correlation 
function X4(fc;£) was monitored in two recent simula- 
tional investigations [H, [24| • Both studies used Newto- 
nian rather than Brownian dynamics and thus our find- 
ings cannot be directly compared to their results. How- 
ever, it is interesting to note that Ref. [13] found that 
the time at which Xi(k',t) reaches its maximum value 
increases as the wave vector decreases from the peak po- 
sition of the structure factor, in rough agreement with 



Fig. [2] The authors of Ref. [24| did not comment on 
the dependence of the maximum value of Xi(k; t) on the 
wave vector but from their Fig. 1 it seems that it does 
not increase with decreasing wave vector, in contrast to 
our Figs. [2] and [3 The second study, Ref. [23|, shows the 
dependence of the maximum value of Xi{k] t) [25j on the 
wave vector. It exhibits a broad maximum located at a 
wave vector slightly smaller than the peak position of the 
static structure factor. Qualitatively, this disagrees with 
our results shown in Figs. [2]and[5l At present, the origin 
of the differences between our findings and the results of 
Refs. [H H|] is unclear. 

In closing, we would like to emphasize that various 
four-point functions quantifying dynamic heterogeneities 
have been studied in computer simulations for more than 
a decade. On the other hand, theoretical predictions for 
these functions started to appear only in the last few 
years. Quantitative comparison between simulations and 
theories is still in infancy. We hope that this work will 
stimulate further effort in this direction. 
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